To do

Organization

Split into three parts:

Projections

Packages & functions

library(plyr)
library(tidyverse)
library(popdemo)
library(popbio)
library(Rcompadre)
library(Rage)

Load functions.

source('R/pop01_param_poun.R')
source('R/pop03_doproj.R')
source('R/pop03_doproj_stoch.R')

Projection parameters.

Each row holds a unique set of parameters for a matrix population model. This reflects diverse recovery (conservation action) and extinction (threat) scenarios.

dt <- pop01_param_poun()
nf <- 10
dt$adultF_n <- nf
ceiling_threshold <- nf + (nf * 0.2)
The first few rows should look like this….
species type first_year a1 a2 a3 a4 b1 b2 b3 b4 c1 c2 c3 c4 d1 d2 d3 d4 akey adultF_n
Podocnemis unifilis base 0.0 0 0 0 10.74 0.0 0.33 0 0 0 0.17 0.29 0 0 0 0.11 0.93 poun_base_0 NA
Podocnemis unifilis base 0.1 0 0 0 10.74 0.1 0.33 0 0 0 0.17 0.29 0 0 0 0.11 0.93 poun_base_0.1 NA
Podocnemis unifilis base 0.2 0 0 0 10.74 0.2 0.33 0 0 0 0.17 0.29 0 0 0 0.11 0.93 poun_base_0.2 NA

Each row can be easily converted for use in a stage-based matrix population model. The following code shows how the first row becomes a stage-based matrix.

stage_names <- c("a1", "a2", "a3", "a4",
                   "b1", "b2", "b3", "b4",
                   "c1", "c2", "c3", "c4",
                   "d1", "d2", "d3", "d4")
vpop <- unlist(dt[3 , stage_names])
pop_mat <- matrix(vpop, byrow = TRUE, ncol=4)
dimnames(pop_mat) <- list(c("a", "b", "c", "d"),
                            c( "a", "b", "c", "d"))
The matrix should look like this:
a b c d
a 0.000000 0.000000 0.000000 10.741500
b 0.200000 0.333333 0.000000 0.000000
c 0.000000 0.166667 0.285714 0.000000
d 0.000000 0.000000 0.114286 0.930000

This matrix includes information of population growth and survival and reproduction.

See https://cran.r-project.org/web/packages/Rage/vignettes/a01_GettingStarted.html

We can show the different matrix components converting to the structure recommended by …

meta <- data.frame(idNum = 3, 
                   SpeciesAccepted = dt[3 , 'species'],
                   type = dt[3 , 'type'], 
                   first_year = dt[3 , 'first_year']) 
stageInfo <- list(
  data.frame(
    MatrixClassOrganized = rep("active", 4),
    MatrixClassAuthor = c("eggs/hatchling", "juvenile-early", "juvenile-late", "adult")
  ))
# Simple full A matrix
l_pop_mat <- list(mpm = pop_mat)
x <- cdb_build_cdb(mat_a = l_pop_mat, metadata = meta, stages = stageInfo)

# Seperate into U, F, C
mat_u1 <- rbind(
  c(dt[3, 'a1'], dt[3, 'a2'], dt[3, 'a3'], 0),
  c(dt[3, 'b1'], dt[3, 'b2'], dt[3, 'b3'], dt[3, 'b4']),
  c(dt[3, 'c1'], dt[3, 'c2'], dt[3, 'c3'], dt[3, 'c4']),
  c(dt[3, 'd1'], dt[3, 'd2'], dt[3, 'd3'], dt[3, 'd4'])
)
mat_f1 <- rbind(
  c(0.0, 0.0, 0.0, dt[3, 'a4']),
  c(0.0, 0.0, 0.0, 0.0),
  c(0.0, 0.0, 0.0, 0.0), 
  c(0.0, 0.0, 0.0, 0.0)
)
mat_c1 <- rbind(
  c(0.0, 0.0, 0.0, 0.0),
  c(0.0, 0.0, 0.0, 0.0),
  c(0.0, 0.0, 0.0, 0.0), 
  c(0.0, 0.0, 0.0, 0.0)
)
l_u <- list(m_u = mat_u1)
l_f <- list(m_f = mat_f1)
l_c <- list(m_c = mat_c1)

my_comadre <- cdb_build_cdb(
  mat_u = l_u, mat_f = l_f, mat_c = l_c,
  metadata = meta, 
  stages = stageInfo
)

matA(my_comadre)
## [[1]]
##      [,1]     [,2]     [,3]    [,4]
## [1,]  0.0 0.000000 0.000000 10.7415
## [2,]  0.2 0.333333 0.000000  0.0000
## [3,]  0.0 0.166667 0.285714  0.0000
## [4,]  0.0 0.000000 0.114286  0.9300
matA(x)
## [[1]]
##     a        b        c       d
## a 0.0 0.000000 0.000000 10.7415
## b 0.2 0.333333 0.000000  0.0000
## c 0.0 0.166667 0.285714  0.0000
## d 0.0 0.000000 0.114286  0.9300
matU(my_comadre)
## [[1]]
##      [,1]     [,2]     [,3] [,4]
## [1,]  0.0 0.000000 0.000000 0.00
## [2,]  0.2 0.333333 0.000000 0.00
## [3,]  0.0 0.166667 0.285714 0.00
## [4,]  0.0 0.000000 0.114286 0.93
matF(my_comadre)
## [[1]]
##      [,1] [,2] [,3]    [,4]
## [1,]    0    0    0 10.7415
## [2,]    0    0    0  0.0000
## [3,]    0    0    0  0.0000
## [4,]    0    0    0  0.0000
matC(my_comadre)
## [[1]]
##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    0    0    0    0
## [3,]    0    0    0    0
## [4,]    0    0    0    0

We can get life history metrics

# the time required for a population to 
# increase by a factor of R0 (the net reproductive rate)
gen_time(matU = matU(my_comadre)[[1]], matF(my_comadre)[[1]]) # 17.15
## [1] 17.15455
# the average parent-offspring age difference # 16.22
gen_time(matU = matU(my_comadre)[[1]], matF(my_comadre)[[1]], method = "age_diff")
## [1] 16.22243
# expected age at reproduction for a cohort 5.69
Rage::gen_time(matU = matU(my_comadre)[[1]], matF(my_comadre)[[1]], method = "cohort")
## [1] 5.696541
Rage::life_expect_mean(matU = matU(my_comadre)[[1]], start = 1)
## [1] 1.484286
sum(Rage::life_expect_mean(matU(my_comadre)[[1]], start = NULL))
## [1] 21.87715

Plot the life cycle diagram.

plot_life_cycle(matA = matA(my_comadre)[[1]], 
                stages = c("eggs/hatchling", "juvenile-early", "juvenile-late", "adult"))

Eggs to adult

Using survival values from base model.

# The probability of reaching reproductive maturity before death for an egg. 
# Expressed as proportion (0 - 1)
# # 0.008000029 # 0.8%
mp_poun <- mature_prob(matU = matU(my_comadre)[[1]], matR = matF(my_comadre)[[1]])

# To produce one reproductive adult need 125 eggs.
1/mp_poun # 125
## [1] 124.9996
# 7 years for a female to lay sufficient eggs 
ceiling(125 / 20)
## [1] 7
t0 <- 1000 # number of eggs
n_t1 <- t0 * 0.2       # 200
n_t2 <- n_t1 * 0.333333 # 66.6
n_t3 <- n_t2 * 0.333333 # 22.2
n_t4 <- n_t3 * 0.285714 #  6.3
n_t5 <- n_t4 * 0.285714 #  1.8
n_t6 <- n_t5 * 0.93     #  1.6
n_t7 <- n_t6 * 0.93     #  1.5
n_t8 <- n_t7 * 0.93     #  1.4
n_t9 <- n_t8 * 0.93     #  1.3
ceiling(t0 / n_t5) # 552 eggs to produce one adult (maturity 5 years)
## [1] 552
ceiling(t0 / n_t6) # 593 eggs to produce one adult (maturity 6 years)
## [1] 593
ceiling(t0 / n_t7) # 638 eggs to produce one adult (maturity 7 years)
## [1] 638
ceiling(t0 / n_t8) # 686 eggs to produce one adult (maturity 8 years)
## [1] 686
ceiling(t0 / n_t9) # 737 eggs to produce one adult (maturity 9 years)
## [1] 737
ceiling(552 / 20) # 28 years for a female to lay sufficient eggs (maturity 5 years)
## [1] 28
ceiling(737 / 20) # 37 years for a female to lay sufficient eggs (maturity 9 years)
## [1] 37

Projection models.

Deterministic models

# project
dout <- plyr::ddply(dt, 
                    c("species", "type", "first_year","akey"), .fun =  pop03_doproj)
dout$arun <- 1

# Model summaries
model_sum <- dout |> 
  group_by(species, type, first_year, lambda, 
           gen_time, gen_age_diff, life_exp, life_exp_adult, mat_prob, eggs_to_adult) |> 
  summarise(fem_t0 = max(fem_t0), 
            fem_min = min(fem),
              fem_max = max(fem)) |> 
  ungroup()
lambda_n <- length(unique(model_sum$lambda)) # 50
lambda_mean <- mean(model_sum$lambda) # 0.9432
lambda_sd <- sd(model_sum$lambda) # 0.1506
lambda_min <- min(model_sum$lambda) # 0.4659
lambda_max <- max(model_sum$lambda) # 1.1539

# Export for future use
saveRDS(dout, "inst/other/dout.rds")

Stochastic models

# Stochastic
#data frame with runs for processing
#nruns <- 100 # 100 gives same pattern as 50
nruns <- 50
dt_stoch <- dt[rep(seq_len(nrow(dt)), nruns), ]
dt_stoch$arun <- rep(1:nruns, each = nrow(dt))
# Approx 90 minutes. 1,212,000 rows. Projections quick. Summaries slow.
# 11:32 - 13:14
dout_stoch <- plyr::ddply(dt_stoch, 
                    c("arun", "species", "type", "first_year","akey"), 
                    .fun =  pop03_doproj_stoch)
table(dout_stoch$model)
table(dout_stoch$type)
model_sum_stoch <- dout_stoch |> 
  group_by(species, type, model, first_year, lambda, lambda_q75,
           gen_time, gen_age_diff_med, gen_age_q75, 
           life_exp_med, life_exp_adult_med, mat_prob_med, mat_prob_q75, 
           eggs_to_adult_med, eggs_to_adult_q75) |> 
  summarise(acount = n(), 
            fem_t0 = max(fem_t0), 
            fem_min = min(fem),
              fem_max = max(fem)) |> 
  ungroup()
# Export for future use
saveRDS(dout_stoch, "inst/other/dout_stoch.rds")

Combine results.

dout <- readRDS("inst/other/dout.rds")
dout_stoch <- readRDS("inst/other/dout_stoch.rds")
# Combine data for plotting
dout_all <- dplyr::bind_rows(dout |> dplyr::select(arun, model, type, first_year, 
                                                   akey, ayear, 
                                                   lambda, gen_time, gen_age_diff, 
                                                   life_exp, life_exp_adult, 
                                                   mat_prob, eggs_to_adult,
                                            fem, fem_t0, fem_diff, change50_flag, 
                                            change30_flag,
                                            double_flag) |> 
                               dplyr::mutate(lambda_lcl = NA, lambda_ucl = NA, 
                                             lambda_sd = NA,  gen_sd = NA), 
                             dout_stoch |> dplyr::select(arun, model, type, first_year, 
                                                         akey, ayear, 
                                                         lambda, lambda_lcl, lambda_ucl, 
                                                         lambda_sd, gen_time, gen_sd, 
                                                         gen_age_diff, 
                                                   life_exp, life_exp_adult, 
                                                   mat_prob, eggs_to_adult,
                                            fem, fem_t0, fem_diff, change50_flag, 
                                            change30_flag,
                                            double_flag))
# Limit adult female number to maximum (20% above original for baseline).
summary(dout_all$fem)
#    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#       0        0        1     1389       11 16560701
dout_all[which(dout_all$fem > ceiling_threshold), 'fem' ] <- ceiling_threshold
# summary(dout_all$fem)
#     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
# 0.000000  0.005763  1.485889  4.522590 11.240127 12.000000 
 
# Factors in right order
dout_all$modelf <- 1
dout_all[which(dout_all$model=="Stochastic uniform") , 'modelf'] <- 2
dout_all[which(dout_all$model=="Stochastic equal") , 'modelf'] <- 3
dout_all[which(dout_all$model=="Stochastic bad x2") , 'modelf'] <- 4
dout_all[which(dout_all$model=="Stochastic bad x4") , 'modelf'] <- 5
dout_all$modelf <- as.factor(dout_all$modelf)
levels(dout_all$modelf) <- c("Deterministic", "Stochastic uniform", 
                          "Stochastic equal", "Stochastic bad x2", 
                          "Stochastic bad x4")
unique(dout_all$modelf)
table(dout_all$modelf)
dout_all$typef <- 1
dout_all[which(dout_all$type=="female-hunt 2.5%") , 'typef'] <- 2
dout_all[which(dout_all$type=="female-hunt 5%") , 'typef'] <- 3
dout_all[which(dout_all$type=="female-hunt 10%") , 'typef'] <- 4
dout_all[which(dout_all$type=="female-hunt 25%") , 'typef'] <- 5
dout_all[which(dout_all$type=="female-hunt 50%") , 'typef'] <- 6
dout_all$typef <- as.factor(dout_all$typef)
levels(dout_all$typef) <- c("base", "female-hunt 2.5%", 
                            "female-hunt 5%",
                          "female-hunt 10%", "female-hunt 25%", 
                          "female-hunt 50%")
table(dout_all$typef)
# first year survival
dout_all$first_yearf <- as.factor(dout_all$first_year)
fylev <- paste("first-year\nsurvival\n", seq(0, 0.9, by = 0.1), sep = "")
levels(dout_all$first_yearf) <- fylev
# Export for future use. 17/6/2024 - 1,218,060 rows 21 columns.
saveRDS(dout_all, "inst/other/dout_all.rds")

Model lookup

Model summaries to link with scenarios. Summaries across stochastic runs for each model.

  • RedList guidelines p38
    For example, upper and lower quartiles of the projected magnitude of the future reduction (i.e., reductions with 25% and 75% probability) may be considered to represent a plausible range of projected reduction

  • Green Status Presence
    Present when population is at least 1 as measured by lower quartile population.

  • Green Status Viability
    Viable when not declining: lower quartile of lambda >= 1. A species is considered viable in a spatial unit if application of the Regional Red List Guidelines to the population in the spatial unit would result in a categorization of ‘Least Concern’ OR ‘Near Threatened and Not Declining’.

  • Green Status Functional
    Functional when population lower quartile is 10 times base level.

# load data. 2 August 2024 with 1218060 rows and 21 columns
dout_all <- readRDS("inst/other/dout_all.rds")
# 101 years, 5 model types (deterministic with 6060, stochastic with 6060 * 50), 
# 6 harvest levels, 10 first year levels
# Make unique model ID. boot mean is same as mean (at least to 6 decimal places)
# 300 projection models.
model_ref <- dout_all |>
  group_by(akey, modelf, typef, first_year, arun, lambda, gen_time, 
           gen_age_diff, life_exp, life_exp_adult, eggs_to_adult) |> 
  summarise(yc = length(unique(ayear))) |> 
  ungroup() |> 
  group_by(akey, modelf, typef, first_year) |>
  summarise(count_runs = length(unique(arun)), 
            count_years = min(yc),
            lambda_mean = mean(lambda), 
            lambda_min = min(lambda), 
            lambda_max = max(lambda), 
            lambda_sd = sd(lambda),
            lambda_boot_lcl = Hmisc::smean.cl.boot(lambda)["Lower"],
            lambda_boot_ucl = Hmisc::smean.cl.boot(lambda)["Upper"], 
            lambda_q25 = quantile(lambda, probs = 0.25, na.rm = TRUE),
            lambda_q75 = quantile(lambda, probs = 0.75, na.rm = TRUE), 
            gen_mean = mean(gen_time), 
            gen_med = median(gen_time),
            gen_min = min(gen_time), 
            gen_max = max(gen_time), 
            gen_sd = sd(gen_time),
            gen_boot_lcl = Hmisc::smean.cl.boot(gen_time)["Lower"],
            gen_boot_ucl = Hmisc::smean.cl.boot(gen_time)["Upper"], 
            gen_q25 = quantile(gen_time, probs = 0.25, na.rm = TRUE),
            gen_q75 = quantile(gen_time, probs = 0.75, na.rm = TRUE), 
            gen_age_mean = mean(gen_age_diff), 
            gen_age_med = median(gen_age_diff),
            gen_age_min = min(gen_age_diff), 
            gen_age_max = max(gen_age_diff), 
            gen_age_sd = sd(gen_age_diff),
            gen_age_boot_lcl = Hmisc::smean.cl.boot(gen_age_diff)["Lower"],
            gen_age_boot_ucl = Hmisc::smean.cl.boot(gen_age_diff)["Upper"], 
            gen_age_q25 = quantile(gen_age_diff, probs = 0.25, na.rm = TRUE),
            gen_age_q75 = quantile(gen_age_diff, probs = 0.75, na.rm = TRUE), 
            life_exp_mean = mean(life_exp), 
            life_exp_med = median(life_exp),
            life_exp_min = min(life_exp), 
            life_exp_max = max(life_exp), 
            life_exp_sd = sd(life_exp),
            life_exp_boot_lcl = Hmisc::smean.cl.boot(life_exp)["Lower"],
            life_exp_boot_ucl = Hmisc::smean.cl.boot(life_exp)["Upper"], 
            life_exp_q25 = quantile(life_exp, probs = 0.25, na.rm = TRUE),
            life_exp_q75 = quantile(life_exp, probs = 0.75, na.rm = TRUE), 
            life_exp_adult_mean = mean(life_exp_adult), 
            life_exp_adult_med = median(life_exp_adult),
            life_exp_adult_min = min(life_exp_adult), 
            life_exp_adult_max = max(life_exp_adult), 
            life_exp_adult_sd = sd(life_exp_adult),
            life_exp_adult_boot_lcl = Hmisc::smean.cl.boot(life_exp_adult)["Lower"],
            life_exp_adult_boot_ucl = Hmisc::smean.cl.boot(life_exp_adult)["Upper"], 
            life_exp_adult_q25 = quantile(life_exp_adult, probs = 0.25, na.rm = TRUE),
            life_exp_adult_q75 = quantile(life_exp_adult, probs = 0.75, na.rm = TRUE),
            eggs_to_adult_mean = mean(eggs_to_adult), 
            eggs_to_adult_med = median(eggs_to_adult),
            eggs_to_adult_min = min(eggs_to_adult), 
            eggs_to_adult_max = max(eggs_to_adult), 
            eggs_to_adult_sd = sd(eggs_to_adult),
            eggs_to_adult_boot_lcl = Hmisc::smean.cl.boot(eggs_to_adult)["Lower"],
            eggs_to_adult_boot_ucl = Hmisc::smean.cl.boot(eggs_to_adult)["Upper"], 
            eggs_to_adult_q25 = quantile(eggs_to_adult, probs = 0.25, na.rm = TRUE),
            eggs_to_adult_q75 = quantile(eggs_to_adult, probs = 0.75, na.rm = TRUE)
            
            ) |> 
  ungroup() |> 
  arrange(typef, first_year, modelf) |>
  mutate(modelid = paste(akey, as.numeric(modelf), sep = "_"), 
         modelkey = row_number()) |> 
  relocate(modelkey, modelid)

# Estimates of generation time for summary.
# Bienvenu and Legendre (2015) https://doi.org/10.1086%2F681104 - 
# defined as the mean age of mothers at birth.
# IUCN definition - is the mean age at which a cohort of individuals produce offspring.
gen_mean <- model_ref |> 
  filter(lambda_q75 >= 1, first_year < 0.5) |> pull(gen_age_mean) |> mean()
gen_median <- model_ref |> 
  filter(lambda_q75 >= 1, first_year < 0.5) |> pull(gen_age_mean) |> median()
gen_q25 <- model_ref |> 
  filter(lambda_q75 >= 1, first_year < 0.5) |> pull(gen_age_mean) |> 
  quantile(probs = 0.25)
gen_q75 <- model_ref |> 
  filter(lambda_q75 >= 1, first_year < 0.5) |> pull(gen_age_mean) |> 
  quantile(probs = 0.75)
gen_3 <- gen_mean * 3
gen_3_q25 <- gen_q25 * 3
gen_3_q75 <- gen_q75 * 3

# Add population changes.
model_ref2 <- model_ref |> left_join(dout_all |>
  filter(ayear %in% c(ceiling(gen_3), ceiling(gen_3_q25), ceiling(gen_3_q75), 100)) |> 
  group_by(akey, modelf, typef, first_year, ayear, arun) |> 
  summarise(fem_t0 = mean(fem_t0), 
            femnew = mean(fem),
            femnew_diff = mean(fem_diff)
  ) |> 
    ungroup() |> 
  group_by(akey, modelf, typef, first_year, ayear) |> 
  summarise(fem_t0 = mean(fem_t0), 
            fem = mean(femnew),
            fem_sd = sd(femnew, na.rm = TRUE), 
            fem_q25 = quantile(femnew, probs = 0.25, na.rm = TRUE),
            fem_q75 = quantile(femnew, probs = 0.75, na.rm = TRUE),
            fem_diff = mean(femnew_diff), 
            fem_diff_sd = sd(femnew_diff, na.rm = TRUE), 
            fem_diff_q25 = quantile(femnew_diff, probs = 0.25, na.rm = TRUE),
            fem_diff_q75 = quantile(femnew_diff, probs = 0.75, na.rm = TRUE)) |> 
  ungroup() |> 
    mutate(ayear = paste("t", ayear, sep="")) |>
  pivot_wider(names_from = ayear, 
               values_from = c(fem, fem_sd, fem_q25, fem_q75, fem_diff, 
                               fem_diff_sd, fem_diff_q25, fem_diff_q75)) |> 
  arrange(typef, first_year, modelf)
)

# Add time to change by 50% and 30% 
model_ref3 <- model_ref2 |> 
  left_join(dout_all |> 
  filter(change50_flag == 1) |>
 group_by(akey, modelf, typef, first_year, arun, 
          change50_flag) |>
  summarise(acount = n(), 
            change50_yf = min(ayear), 
            change50_yl = max(ayear)) |> 
  ungroup() |> 
  group_by(akey, modelf, typef, first_year) |> 
  summarise(change50_ymean = floor(mean(change50_yf)), 
            change50_ysd = sd(change50_yf), 
            change50_yq25 = floor(quantile(change50_yf, probs = 0.25, na.rm = TRUE)),
            change50_yq75 = floor(quantile(change50_yf, probs = 0.75, na.rm = TRUE))) |> 
   ungroup()
) |> left_join(dout_all |> 
  filter(change30_flag == 1) |>
 group_by(akey, modelf, typef, first_year, arun, 
          change30_flag) |>
  summarise(acount = n(), 
            change30_yf = min(ayear), 
            change30_yl = max(ayear)) |> 
  ungroup() |> 
  group_by(akey, modelf, typef, first_year) |> 
  summarise(change30_ymean = floor(mean(change30_yf)), 
            change30_ysd = sd(change30_yf), 
            change30_yq25 = floor(quantile(change30_yf, probs = 0.25, na.rm = TRUE)),
            change30_yq75 = floor(quantile(change30_yf, probs = 0.75, na.rm = TRUE))) |> 
   ungroup()
)

# Export for future use
saveRDS(model_ref3, "inst/other/model_lookup.rds") 
write.csv2(model_ref3, "inst/other/model_lookup.csv", row.names = FALSE) 
# summaries of  10050 runs.
#model_all_sum <- dout_all |> 
#  group_by(akey, modelf, typef, first_year, arun) |> 
#  summarise(count_years = length(ayear), 
#            count_years_unique = length(unique(ayear))) |> 
#  ungroup()

Compare lambda. Use this as the basis to establish number of runs to use.

model_lookup <- readRDS("inst/other/model_lookup.rds")

model_lookup |> 
  mutate(recover_flag = if_else(fem_diff_t100 > 0, 1, 0)) |> 
  group_by(recover_flag) |> 
  summarise(model_count = n(),
            lambda_avg = mean(lambda_mean), 
            lambda_min = min(lambda_mean), 
            lambda_max = max(lambda_mean))
## # A tibble: 2 × 5
##   recover_flag model_count lambda_avg lambda_min lambda_max
##          <dbl>       <int>      <dbl>      <dbl>      <dbl>
## 1            0         212      0.842      0.466      0.995
## 2            1          88      1.04       0.973      1.15
# 18 where mean lambda is less than one but population increases.
# 100, 144 and 198 with lowest lambda. Results uncertain
model_lookup |> 
  filter(fem_diff_t100 > 0, lambda_mean < 1) |> 
  arrange(lambda_mean)
## # A tibble: 19 × 102
##    modelkey modelid         akey  modelf typef first_year count_runs count_years
##       <int> <chr>           <chr> <fct>  <fct>      <dbl>      <int>       <int>
##  1      193 poun_female-hu… poun… Stoch… fema…        0.8         50         101
##  2      144 poun_female-hu… poun… Stoch… fema…        0.8         50         101
##  3      100 poun_female-hu… poun… Stoch… fema…        0.9         50         101
##  4      198 poun_female-hu… poun… Stoch… fema…        0.9         50         101
##  5      149 poun_female-hu… poun… Stoch… fema…        0.9         50         101
##  6       35 poun_base_0.6_5 poun… Stoch… base         0.6         50         101
##  7       84 poun_female-hu… poun… Stoch… fema…        0.6         50         101
##  8       73 poun_female-hu… poun… Stoch… fema…        0.4         50         101
##  9       40 poun_base_0.7_5 poun… Stoch… base         0.7         50         101
## 10      133 poun_female-hu… poun… Stoch… fema…        0.6         50         101
## 11       24 poun_base_0.4_4 poun… Stoch… base         0.4         50         101
## 12       89 poun_female-hu… poun… Stoch… fema…        0.7         50         101
## 13       45 poun_base_0.8_5 poun… Stoch… base         0.8         50         101
## 14       50 poun_base_0.9_5 poun… Stoch… base         0.9         50         101
## 15       94 poun_female-hu… poun… Stoch… fema…        0.8         50         101
## 16       18 poun_base_0.3_3 poun… Stoch… base         0.3         50         101
## 17       29 poun_base_0.5_4 poun… Stoch… base         0.5         50         101
## 18       78 poun_female-hu… poun… Stoch… fema…        0.5         50         101
## 19      138 poun_female-hu… poun… Stoch… fema…        0.7         50         101
## # ℹ 94 more variables: lambda_mean <dbl>, lambda_min <dbl>, lambda_max <dbl>,
## #   lambda_sd <dbl>, lambda_boot_lcl <dbl>, lambda_boot_ucl <dbl>,
## #   lambda_q25 <dbl>, lambda_q75 <dbl>, gen_mean <dbl>, gen_med <dbl>,
## #   gen_min <dbl>, gen_max <dbl>, gen_sd <dbl>, gen_boot_lcl <dbl>,
## #   gen_boot_ucl <dbl>, gen_q25 <dbl>, gen_q75 <dbl>, gen_age_mean <dbl>,
## #   gen_age_med <dbl>, gen_age_min <dbl>, gen_age_max <dbl>, gen_age_sd <dbl>,
## #   gen_age_boot_lcl <dbl>, gen_age_boot_ucl <dbl>, gen_age_q25 <dbl>, …
# 10 where t100 females can be more or less than t0.
model_lookup |> 
  mutate(flag_uncertain = if_else((fem_q25_t100 < fem_t0) & (fem_q75_t100 > fem_t0), 
                                  1, 0)) |> 
  filter(flag_uncertain == 1) |> 
  arrange(lambda_mean)
## # A tibble: 9 × 103
##   modelkey modelid          akey  modelf typef first_year count_runs count_years
##      <int> <chr>            <chr> <fct>  <fct>      <dbl>      <int>       <int>
## 1      193 poun_female-hun… poun… Stoch… fema…        0.8         50         101
## 2      144 poun_female-hun… poun… Stoch… fema…        0.8         50         101
## 3      100 poun_female-hun… poun… Stoch… fema…        0.9         50         101
## 4      128 poun_female-hun… poun… Stoch… fema…        0.5         50         101
## 5       35 poun_base_0.6_5  poun… Stoch… base         0.6         50         101
## 6       84 poun_female-hun… poun… Stoch… fema…        0.6         50         101
## 7       73 poun_female-hun… poun… Stoch… fema…        0.4         50         101
## 8       40 poun_base_0.7_5  poun… Stoch… base         0.7         50         101
## 9       24 poun_base_0.4_4  poun… Stoch… base         0.4         50         101
## # ℹ 95 more variables: lambda_mean <dbl>, lambda_min <dbl>, lambda_max <dbl>,
## #   lambda_sd <dbl>, lambda_boot_lcl <dbl>, lambda_boot_ucl <dbl>,
## #   lambda_q25 <dbl>, lambda_q75 <dbl>, gen_mean <dbl>, gen_med <dbl>,
## #   gen_min <dbl>, gen_max <dbl>, gen_sd <dbl>, gen_boot_lcl <dbl>,
## #   gen_boot_ucl <dbl>, gen_q25 <dbl>, gen_q75 <dbl>, gen_age_mean <dbl>,
## #   gen_age_med <dbl>, gen_age_min <dbl>, gen_age_max <dbl>, gen_age_sd <dbl>,
## #   gen_age_boot_lcl <dbl>, gen_age_boot_ucl <dbl>, gen_age_q25 <dbl>, …

Plots

Compare lambda.

mean_values <- model_ref |> 
  group_by(typef) |> 
  summarise(l_mean = mean(lambda_mean), 
            l_lcl = Hmisc::smean.cl.boot(lambda_mean)["Lower"], 
            l_ucl = Hmisc::smean.cl.boot(lambda_mean)["Upper"], 
            )

fig_model_lambda <- model_ref |> 
  ggplot(aes(x = first_year, y = lambda_mean, colour = modelf)) + 
  geom_hline(data = mean_values, aes(yintercept = l_mean), 
             colour = "black", linewidth = 0.6,
             linetype = "dashed") + 
    geom_hline(data = mean_values, aes(yintercept = l_lcl), 
             colour = "grey", linewidth = 0.6,
             linetype = "dashed") + 
      geom_hline(data = mean_values, aes(yintercept = l_ucl), 
             colour = "grey", linewidth = 0.6,
             linetype = "dashed") + 
  geom_hline(yintercept = 1, 
             colour = "magenta", linewidth = 0.4) + 
  geom_point() + 
  stat_smooth(se = FALSE) + 
  scale_x_continuous("first year survival and graduation", 
                     breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_color_viridis_d() +
  facet_wrap(~ typef, ncol = 6) + 
  theme(legend.position="top") + 
  guides(colour=guide_legend(nrow=2,byrow=TRUE)) + 
  labs(y = "lambda", colour = "model")
# save
png(file = "inst/other/fig_unifilis_lambda.png", bg = "transparent", 
    type = c("cairo"), 
    width = 9, height = 4, units = "in", res=600)
fig_model_lambda
invisible(dev.off())
Check plot.
Model lambda for *Podocnemis unifilis*

Model lambda for Podocnemis unifilis

Compare number of adult females

fig_proj <- dout_all |> 
  ggplot(aes(x = ayear, y = fem, colour = modelf)) + 
  geom_point(size=0.1, alpha=0.2) + 
  stat_smooth(se = FALSE) + 
  scale_colour_viridis_d("model") + 
  scale_y_continuous(limits = c(0, ceiling_threshold), 
                     breaks =  seq(2, 12, by = 2)) +
  labs(y="Reproductive females",
       x="Time (years)",
       ) +
  theme_bw() +
  theme(plot.title.position = "plot") + 
  facet_grid(first_yearf ~ typef)

# save
png(file = "inst/other/fig_unifilis_projections.png", bg = "transparent", 
    type = c("cairo"), 
    width = 10, height = 8, units = "in", res=600)
fig_proj
invisible(dev.off())
Check plot.
Model scenarios for *Podocnemis unifilis*

Model scenarios for Podocnemis unifilis

Number of adult females along 80 km2 of river.

# Here use density values from "tracaja_dist_5km_4z_beforeafter.R" to establish "impact"
# Adult population before - after
an_b4 <- ceiling((80 * 1.0035150)) # before = 81 in 80 km2 of river
an_b4_lci <- ceiling((80 * 0.38838535))
an_b4_uci <- ceiling((80 * 2.5928949))
an_aft <- ceiling((80 *  0.1542282)) # after  = 13 in 80 km2 of river
an_aft_lci <- ceiling((80 * 0.04021797))
an_aft_uci <- ceiling((80 * 0.5914359))